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Abstract 

We present a perfect simulation algorithm for measures that are absolutely continu- 
ous with respect to some Poisson process and can be obtained as invariant measures 
of birth-and-death processes. Examples include area- and perimeter- interacting point 
processes (with stochastic grains), invariant measures of loss networks, and the Ising 
contour and random cluster models. The algorithm does not involve couplings of 
the process with different initial conditions and it is not tied up to monotonicity 
requirements. Furthermore, it directly provides perfect samples of finite windows of 
the infinite-volume measure, subjected to time and space "user-impatience bias". 
The algorithm is based on a two-step procedure: (i) a perfect-simulation scheme for 
a (finite and random) relevant portion of a (space-time) marked Poisson processes 
(free birth-and-death process, free loss networks), and (ii) a "cleaning" algorithm 
that trims out this process according to the interaction rules of the target process. 
The first step involves the perfect generation of "ancestors" of a given object, that 
is of predecessors that may have an infiuence on the birth-rate under the target 
process. The second step, and hence the whole procedure, is feasible if these "ances- 
tors" form a finite set with probability one. We present a sufficiency criteria for this 
condition, based on the absence of infinite clusters for an associated (backwards) 
oriented percolation model. The criteria is expressed in terms of the subcritical- 
ity of a majorizing multi-type branching process, whose corresponding parameter 
yields bounds for errors due to space-time "user-impatience bias" . The approach has 
previously been used, as an alternative to cluster expansion techniques, to extract 
properties of the invariant measures involved. 
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1 Introduction 



Perfect simulations or exact sampling are labels for a recently developed set 
of techniques designed to produce output whose distribution is guaranteed 
to follow a given probability law. These techniques are particularly useful in 
relation with Markov Chain Monte Carlo, and their range of applicability is 
rapidly growing (see pittp : //dimacs . rutgers . edu/~dbwilson/exact| ). 



There are several techniques for perfect simulation of Markov processes. The 
most popular ones can be classified in two categories: Propp and Wilson's Cou- 
pling from the Past (CFTP) and Fill's Interruptible Algorithm. The first type 
of technique applies, in its original version, to invariant measures of Markov 
processes with a finite coalescence time. That is, of processes for which there 
exists a coupling among trajectories such that with probability one the tra- 
jectories starting from all possible initial states coalesce in a finite time. This 
includes all irreducible Markov processes with a finite state space. The coa- 
lescence property becomes difficult to check if the state space is very large. 
The problem can be overcome for processes with the following monotonicity 
property: there must exist a "maximal" and a "minimal" state and a coupling 
such that the coalescence of coupled trajectories starting from these two states 
implies the coalescence of all other trajectories ("monotone coupling"). Exam- 
ples of processes with this property include Glauber dynamics of spin systems 
with the FKG property (Propp and Wilson, 1996). Other perfect simulation 
techniques are based on backward coupling of embedded regeneration times 
(Corcoran and Tweedie, 2001), finitary coding (Haggstrom and Steif, 2000; 
van den Berg and Steif, 1999), tempering algorithms, cluster representation 
of a Markov chain (Cai, 1999 preprint), regenerative construction (Comets, 
Fernandez and Ferrari, 2002), embedding the problem into the coloring of a 
graph (Fill and Huber, 2000) , representation as an infinite mixture (Hobert 
and Robert, 2000) among others. 



2 



The basic CFTP algorithm, sometimes called vertical CFTP, is in general not 
applicable to processes with infinite state space. Indeed, most of them lack 
uniform ergodicity, a property shown (Foss and Tweedie, 1998) to be equiv- 
alent to the existence of a coalescence scheme as above. To cope with this 
situation, Kendall (1997 and 1998) introduced dominated CFTP (also called 
horizontal CFTP and coupling into and from the past) . This extension also 
requires the state space to have a partial order, as well as the existence of a 
monotone coupling among the target process and two reversible sandwiching 
processes, which must be easy to sample. Algorithms of this type are avail- 
able for attractive point processes and, through a minor modification, also for 
repulsive point processes (Kendall, 1998). Similarly, Haggstrom, van Lieshout 
and M0llcr (1999) combined ideas from CFTP and the two-component Gibbs 
sampler to perfect simulate from process in infinite spaces which do not have 
maximal (or minimal) elements. 

The interruptible algorithm proposed by Fill (1998) (see also Thonnes, 1997) is 
an acceptance-rejection scheme which applies to invariant measures of Markov 
processes whose time-reversed process has a monotonicity property. Thus its 
range of applicability overlaps with that of the CFTP algorithm at reversible 
monotone processes like Glauber dynamics of attractive automata or ferromag- 
netic spin systems and attractive point processes. Later developments have 
made Fill's algorithm applicable to other processes as well (Fill, Machida, 
Murdoch and Rosenthal, 2000). An important advantage of this algorithm is 
that it is free of the so called impatient-user bias: no bias is introduced if the 
user aborts a long run of the algorithm. 

Kendall (1997, 1998) and Kendall and M0ller (2000) proposed dominated 
CFTP schemes applicable to finite- volume measures which are absolutely con- 
tinuous with respect to a finite Poisson point process and that can be obtained 
as the invariant measure of an interacting spatial birth-and-death process. 
These algorithms are based on two ingredients: (i) the "thinning" of a space- 
time marked Poisson process, and (ii) the coupled construction of upper and 
lower processes whose coalescence signals the output of a perfect sample. See 
the recent review of M0ller (2000) for more references. 

In this paper we propose a new perfect-simulation algorithm which applies to 
the same type of measures but has the following distinctive features: 



3 



• We sample directly from a time stationary realization of the process. There 

is no coalescence criterion, either between coupled realizations or between 
sandwiching processes. The scheme neither requires nor takes advantage of 
monotonicity properties. 

• The scheme directly samples a finite window of the equilibrium measure in 
infinite-volume. In contrast, Kendall (1997, 1998) focus on finite windows 
with fixed boundary conditions, and the infinite-volume limit requires an 
additional process of "perfect simulation in space". Our construction is in 
the spirit of the algorithms proposed by van den Berg and Steif (1999) 
and by Haggstrom and Steif (1999) to simulate infinite- volume measures 
for nearest neighbor interactions in a lattice at high temperature or "high 
noise". In a sense, our algorithm is complementary to those, because it 
applies to regimes where they break down (e.g. at low-temperature). We 
point out that before the arrival of the perfect simulation wave, Ferrari 
(1990), van den Berg (1993) and van den Berg and Maes (1994) have also 
proposed construction schemes for (infinite- vohime) Gibbs measures of spin 
systems that can be easily transcribed into perfect-simulation algorithms. 

• The construction has the added value of being a proven theoretical tool for 
the analysis of properties of the target measure. For instance, in Fernandez, 
Ferrari and Garcia (1998 and 2001) we used it to obtain mixing proper- 
ties, finite-volume corrections and the asymptotic (in temperature) distri- 
bution of "defects" of the low-temperature Ising translation-invariant ex- 
tremal measures. Ferrari and Garcia (1998) used a similar construction to 
show ergodicity of a family of loss networks in M. 

• More generally, the construction can be used as an alternative to the cluster- 
expansion technology (Brydges 1986, Kotecky-Preiss 1986, Dobrushin 1996) 
for the study of spin models, at least those with a flipping symmetry. In 
fact, it seems to have a region of validity more extended than usual cluster- 
expansion approaches. 

Our algorithm does involve the "thinning" of a marked Poisson process — the 
free process — which dominates the birth-and-death process, and it involves 
a time-backward and a time-forward sweep. But these procedures are per- 
formed in a form quite different from previous algorithms. The initial stage 
of our construction is done towards the past, starting with a flnite window 
and retrospectively looking to ancestors, namely to those births in the past 
that could have (had) an influence on the current birth. The construction of 
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the clan of ancestors constitutes the time-backward sweep of the algorithm. 
Once this clan is completely constructed, the algorithm proceeds in a time- 
forward fashion "cleaning up" successive generations according to appropriate 
penalization schemes. This "ancestors approach" offers some noteworthy ad- 
vantages: 

(i) The algorithm constructs only the portion of the birth-and-death process 
strictly needed for the final window at t = 0. This economy has two impor- 
tant consequences: First, we can sample directly from the infinite- volume 
measure, without boundary effects. Second, the scheme works for point pro- 
cesses with quite general grain distribution, for instance chosen from an 
unbounded family of objects. As a consequence, it can be applied to the 
simulation of loss networks or of Peierls contours of the low-temperature 
Ising model. 

(ii) Perfect samphng is assured once the algorithm determines the "first" ances- 
tors, that is those ancestors that themselves do not have ancestors. Thus, 
the algorithm determines by itself, in a single sweep, how far back into 
the past the simulation must go. This contrasts with usual CFTP schemes 
where algorithms may have to be iterated several times, going further and 
further into the past, until coalescence is achieved. 

The relation "being ancestor of" induces a backwards in time contact/ oriented 
percolation process. The algorithm is applicable as long as this oriented per- 
colation process is subcritical. This imphes the following hmitations of our 
scheme: 

• It works at low density of objects, at least in infinite volume. It may work at 
higher densities in finite volume, but we have not pursued this investigation. 

• The birth-rate of objects must be uniformly bounded. This is necessary to 
guarantee the existence of the dominating free birth-and-death process. 

Our algorithm does not rely on any type of monotonicity. Therefore for mono- 
tone systems our algorithm probably needs to go further back into the past 
than dominated CFTP or other schemes that exploit monotonicity. This loss 
in efficiency could be compensated, at least partially, by the "economy" and 
"single-sweep" features mentioned in (i) and (ii) above. As Prof. Kendall 
pointed out (private communication), dominated CFTP constructions "are 
wasteful in that they simulate past grains without regard to whether or not 
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they are in the relevant percolation cluster, but efficient in that they use some 
kind of monotonicity to detect whether or not one needs to investigate further 
back into the past." As a counterpart, insensitivity to monotonicity amounts 
to generality and versatility. In particular, this versatility could be used to 
offset the limitation imposed by the low-density constraint. 

As an illustrative analogy, let us present a parallel with what happens in 
statistical mechanics, where studies usually rely on two types of methods: (1) 
those based on correlation inequalities, and (11) those expansion based. (These 
are not the only methods, others include exact solutions and more abstract 
arguments based on compactness or convexity.) Types (I) and (II) are mutu- 
ally complementary. Correlation inequalities yield very strong results, often 
valid over whole regions in parameter space. Nevertheless, these results are 
rigid in that their validity depends on very precise symmetry (monotonicity) 
properties that can be easily destroyed by even infinitesimal perturbations. 
In contrast, expansion-based techniques are very versatile and robust. While 
their a priori range of convergence is limited — it is restricted to low densities 
or high temperatures — , it is often the case that suitable changes of variables 
place other regions of parameter space within the scope of expansion methods. 
For instance, expansion-based studies of high-temperature spin models work 
with (interacting) random walks (see, e.g. Dobrushin, 1996). Closer to, and 
above, the critical temperature alternative expansions are available (Olivieri 
and Picco, 1990; Fernandez, Frohlich and Sokal, 1992). At low temperature 
the right variables are the contours, and to get closer to (but below) the crit- 
ical temperature coarse grained contours are needed (Gaw§dzki, Kotecky and 
Kupiainen, 1987). In fact, the belief is that there always are "good" variables 
that make everything diluted enough. 

Existing dominated-CFTP algorithms are reminiscent of methods of type (I). 
They are very effective and apply for large intervals of rates. But they are also 
very specific, small alterations in the models could affect monotonicity and 
render an algorithm inapplicable. Our scheme could, perhaps, play a role sim- 
ilar to the studies of type (II) for simulation purposes. In this sense, it is crucial 
that monotonicity requirements be absent from the procedure. The goal is to 
change variables so to fall into a low-rate Poissonian (birth-and-death) process. 
Such changes will in general destroy any (obvious) monotonicity property. As 
an example, our scheme is capable to deal with Peierls contours, and hence to 
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provide an exact-sampling algorithm for the low-tempcraturc Ising model (of 
course, it is "exact" modulo time and space user-impatience) . This is a region 
inaccessible to pre-existing algorithms. The use of other random objects (see 
the end of the previous paragraph) could yield analogous algorithms for other 
regions of the phase diagram. 

The comparison of our algorithm with expansion methods is, in fact, more than 
just an analogy. Its theoretical basis has been used to construct an alternative 
to usual expansion methods in statistical mechanical (Fernandez, Ferrari and 
Garcia, 1998 and 2001). This alternative has a provable region of validity that 
exceeds that of usual cluster-expansion treatments. 

For the sake of completeness we start with the definitions of the most con- 
spicuous space processes whose distributions we can perfect-simulate (Section 
2). Examples include area- and perimeter-interacting point process, invariant 
measures of loss networks, the random cluster model and the contour repre- 
sentation of the or "— " Ising measures at low temperature. Its relation 
with birth-and-death processes is discussed immediately after (Section 3), to- 
gether with the basic simulation approach for the latter. The perfect simulation 
scheme is finally presented in Section 5. Its central piece is the time-backward 
construction of the clan of ancestors of a Poissonian birth-and-death process. 



2 Point processes 



Let G be a measurable space and v a Radon measure on G. Typically G is 
M'^, Z'^, R'' X G' or If- x G', where G' is a set of "animals" or "marks" . Let 



5 = e N° : ^(7) > only for a countable set of 7 e G} 



A point process is a random element N ^ S. We denote with the law of a 
point process A^. N is interpreted either as a random configuration of points 
or a random counting measure on G. 
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Poisson Process The first example is a Poisson process on G with intensity 
measure u. Its law is characterized by 

/(AT : N{B) = A;) = e-"(^V(S)VA;! 

for measurable B C G; besides, under //° N{Bi) are independent if Bi are 
disjoint. When G = R*^ and A C we call ix\ the law of n A. We 
call a Poisson process on G = M*^ homogeneous when z/(A) is a function of 
£(A), the Lcbcsgue measure of A. Similarly, when G = Z'^, the process is called 
homogeneous when is a function of |A|, the number of points in Z^'flA. In 
this case, the intensity is proportional to the Lebesgue (respectively, counting) 
measure and the factor of proportionality is called the rate which equals i^(A) 
for any A with unit Lebesgue measure (resp. counting measure). 

Finite total rate. For future purposes we consider the case x M"*") < oo; 
we interpret the last coordinate as time. One can compute the distribution of 
the (not necessarily finite) time ri, the smallest time-coordinate of the points 
(if any) of the process. Indeed, calling A^" the point Poisson process with rate 
I/, for < t < oo, 

P(ti >t) ^ ¥{N{W^ X [0,i)) = 0) = exp(-i/(M'^ x [0,i)) . (2.1) 

In the case of one-dimensional processes {d — 0) the above reads 

P(ti >t) ^ P(A^([0,t)) = 0) = exp(-i/[0,t)) . (2.2) 

In this paper we consider only point processes that are absolutely continuous 
with respect to a Poisson process with law The law of these processes is 
characterized by 

where ^ is the Radon-Nikodim derivative of fj, with respect to 

A Poisson process that appears in the literature is the germ-grain Poisson 
process. In this case G = R*^ x jB°(M'^), where B°(]R'*) is the set of compact 
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Borel sets of R'^. For each alive germ x e M^, g e B^{W^) is the associated grain. 
Assume that the grains are determined by a random variable independent of 
the rest, given by a certain probability distribution tt^, which may depend on 
the germ location x. The intensity u is defined by 

^{d{x,g)) = f{x)'K^{dg)dx. (2.3) 
where f{x) is the intensity of germs. 

Area-interaction point processes These processes have been introduced 
by Baddeley and Van Lieshout (1995). This is a germ-grain process as defined 

above, but the grain shape is fixed and equal to a compact convex G C M*^. 
We only need to keep track of the germs, so G = W^. The intensity u is 
defined by u^dx) = ndx, k is a positive real number. The intersections of the 
grains determine a weight that corrects the otherwise Poissonian distribution 
of germs. The process is absolutely continuous with respect to the Poisson 
process fi^ with intensity u. The law of the area-interaction process is defined 
for bounded windows A C M'' by 

j^-md{N®G) 

I^AidN) = ^ . l^lidN) , (2.4) 
Za[k,(P) 

where /x° is the law of the unit Poisson process in the box A, is a positive 
parameter, ^a(0) is a normalizing constant and A^® G is the coverage process 
given by 

ATeG := [j {x + G}. (2.5) 

xeN 

Strauss process The setup is the same as the area interaction process, but 
now the unit Poisson process is weighted according to an exponential of the 
number of pairs of points closer than a fixed threshold r. The measure is 
defined by 

^^{dN) = ±e'^i^(^)+^^^(^'^) fxlidN) (2.6) 
Za 
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where S{N,A) is the number of unordered pairs such that \\xi — Xj\\ < r. 
The case P2 > ^ was introduced by Strauss (1975) to model the clustering 
of Californian red wood seedling around older stumps. However, (2.6) is not 
integrable in that case (see Kelly and Ripley (1976)). 



Low-temperature Ising model The well-known Peierls contours allow 
to map the "+" or "— " measures of the ferromagnetic Ising model at low 
temperature into an ensemble of objects — the contours — interacting only by 
perimeter-exclusion. See, for instance. Section 5B of Dobrushin, 1996, for a 
concise and rigorous account of this mapping. The (discrete) set G consists 
of contours; these are hypersurfaces formed by a finite number of {d — 1)- 
dimensional unit cubes — links for = 2, plaquettes for higher dimensions — 
centered at points of Z*^ and perpendicular to the edges of the dual lattice Z'* + 
(|, ■ ■ ■ , I). To each contour one can assign an "origin" in Z'^ and say that two 
contours are equivalent if they coincide after a translation of the origin. Calling 
G' the set of contours modulus this class of equivalence, the set G can be 
expressed by G = Z'' x G', where the first coordinate corresponds to the origin 
and the second to the "shape" of the contour. Call two plaquettes adjacent if 
they share a. {d — 2)-dimensional face. A set of plaquettes, 7, is connected if 
for any two plaquettes in 7 there exists a sequence of adjacent plaquettes in 7 
joining them. The set 7 is closed if every (d — 2)-dimensional face is covered 
by an even number of plaquettes in 7. Contours are connected and closed sets 
of plaquettes. For example, in two dimensions contours are closed polygonals. 
Two contours 7 and 9 are said to be compatible if no plaquette of 7 is adjacent 
to a plaquette of 9. In two dimensions, therefore, contours are compatible if 
and only if they do not share the endpoint of a link. In three dimensions two 
compatible contours can share vertices, but not sides of plaquettes. Ising spin 
configurations in a bounded region with "+" (or "— ") boundary condition are 
in one-to-one correspondence with families of pairwise compatible contours. 

Let the compatibility matrix / : G x G — > {0, 1} be defined by 

0, if 7 and 9 are compatible 



(2.7) 
1, otherwise 



10 



The "Poisson process" //^ in is the product of Poisson measures whose 
7-marginal is Poisson with mean 

w(7) := exp(-/?|7|) 

for 7 C A. Here I7I stands for the number of plaquettes of 7. The intensity 
measure v is discrete: dv = f^B) = J2'y£B'>^{l)- 

Let jiA be the measure defined by: for ^ e {0, 1}^, such that ^(7) < 1(7 C A}, 

/^a(0=^( n [1-^(7,^)0/^^(0 (2.8) 
^ 7,e:?(7)?(e)=i 

where /3 is a positive parameter called inverse temperature. The factor Z\ is 
just the normalization. 

The random cluster model Consider A C Z'^ and let M(A) := {(x.y) E 
A X A : \x — y\ = 1} the set of bonds of A. A bond configuration C, G {0, l}"^^) 
is a function from B(A) — > {0, 1}. Bonds assigned 1 are called open, otherwise 
closed. A cluster of C is a set of sites connected with open bonds; sites sur- 
rounded only by closed bonds are clusters of size 1. Let p e [0, 1] and q > 
be parameters and define the finite volume measure 

^A(C) = ^^P°(^)(l-p)^«)g^(^) (2.9) 
^a[P, Q) 

where 0{Q is the number of open bonds of ( in A, C(C) is the number of 
closed bonds and L{() is the number of clusters. The constant Z/^{p, q) is the 
normalization. In other words, ipfy is absolutely continuous with respect to 
the product measure on B(A) with parameter p, with Radon-Nikodim deriva- 
tive q^^^^yz^ip, q). This model has been introduced by Fortuin and Kasteleyn 
(1972); a review can be found in Grimmett (1995). 

Taking the connected sets of bonds as the basic objects, this model can be 
written as in (2.8). More precisely, if one says that two sets of bonds are 
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incompatible whenever they share some vertices, and takes G = {7 C B(A) : 7 
is finite and connected}, the probabihty weights of the model can be written 
in the form (2.8) with the Poisson means 

w{^) = (^) (-) . (2.10) 

Here ^(7) is the number of bonds of 7 and ^^(7) the number of vertices in 
the extremes of the bonds of 7. [That is, 1^(7) = G A such that {x, y) e 
7 or (y.x) E 7 for some y e 7}.] The transformation Y : {0,1}'^*^^^ — > 
{0, 1}^ defined by 

y (C)(7) = 1 <^=^ 7 is a maximally connected set of open bonds of ( 
satisfies 

99a(C) = /XA(y(C)). (2.11) 



3 Birth-and-death processes 

3.1 Definition and examples 

The common feature linking all the spatial processes described in the previous 
section is that all these distributions can be realized as invariant measures of 
spatial interacting birth-and-death processes. 

We consider the state space S of point configurations on G with a Radon 
measure u as in Section 2. 

The free birth death process is characterized by the fact that individuals are 
born at intensity u and last for a random time exponentially distributed of 
mean one. The generator of the free process is 



12 



A'F{n)^lu{d^) [F{rj + S,)-F{r])] 

G 

+ E vh)[F{v-S,)-F{rj)] (3.1) 

7eG:?7(7)>0 

Here is the configuration with only one point at 7 and (77 + ^){0) = f]{0) + 
^{6) (coordinatewise sum). The invariant (and reversible) measure for the free 
process is the Poisson process with intensity u. 

Let be a measure absolutely continuous with respect to with Radon- 
Nikodim derivative ^: dii{r)) — '^{r])d^^{r]). Define 



AF{r,) = J u{d^) [F{V + S,) - F{v)] 



+ E V{1)[F{V-S,)-F{n)] (3.2) 

7eG:»?(7)>0 

The difference with the free process is that in the interacting process the rate 
of birth i'{d'~f) is corrected with the quotient \E'(?7 + S^)/'^{r]), while the rate 
of death remains unaltered. The measure is reversible for the process with 
generator A. To better interpret this dynamics assume 

A$ := sup — — — — < 00 
7,,7 W(?7) 



and define M : G x 5 ^ [0, 1] by 



If 1/ = A^i/, the generator can be rewritten as 
AF{rj) = I P(d7) M(7|77) [F{rj + S^) - F{n)] 

G 

+ E v{l)[F{v-S,)-F{rj)] (3.4) 

7GG:r?(7)>0 
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This dynamics has the following interpretation. When the current configu- 
ration of objects is C,, object 7 attempts to be born with rate u{dj) and is 
effectively born with probability Mi^jl^). The death rate of any object is one. 

The interaction M induces naturally the notion of incompatibility between 
individuals. This is a not necessarily symmetric matrix / : G x G ^ {0, 1} 
defined by 

M(7|O-M(7|e + <^e)|}>0}, (3.5) 



7(7,^) l{sup{ 



where Se is the configuration having unique individual 6 and the supremum 
is taken over the set of those ^ such that ^ and C + ^0 ^-^e in the set of 
configurations (either {0,1}^ or N^). The function I{'y,9) indicates which 
individuals 9 may have an influence in the birth-rate of the individual 7. In 
the case of the spatial point processes described above the matrix is symmetric 
and givenby /((a;,5(), (a;',5(')) = l{{x + g)n{x' + g') ^ 0}. If 1(7, 6*) = 1, — that 
is, if the presence/ absence of 9 modifies the rate of birth of 7 — we say that 9 
is incompatible with 7. For the Ising and random-cluster models one recovers 
(2.7). 



Area interaction Point processes In the repulsive (0 < 1) point process 
(2.4) we have = ^^-"^^(G) and 

u{dx) = Kcp-'^'^^^^dx , (3.6) 

M{x\$,) = 0"^ci(G)-md((x+G)\(^eG)) ^3 

For the attractive (0 > 1) case, 

u{dx) = Kdx (3.8) 

M{x\0 = (p-^Mx+ameG)) _ (3_g) 



Continuous unbounded one- dimensional loss network A loss network 
models, for instance, the occurrence of calls in a communication network. Kelly 
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(1991) reviews several discrete regimes and introduces the following continuous 
generalization. Callers are arranged along an infinitely long cable and each call 
between two points Si, S2 G M on the cable involves just the segment between 
them. The cable has the capacity to carry simultaneously up to C calls past 
any point along its length. Hence, a call attempt between si and S2 G 
Si < S2, is lost if past any point of the interval [si,S2] the cable is already 
carrying C calls. Calls are attempted with initial (leftmost) point following a 
space-time Poisson process with intensity f{x)dx, and (space) lengths given 
by a distribution tt, independent of its leftmost point, with finite mean pi. 
The holding time of a call has exponential distribution with mean one. The 
location of a call, its length and its duration are independent. 

In this case, the germs (x) are the leftmost points of calls and the grains (g) 
are segments with random lengths. This process can be viewed as a spatial 
birth and death process where we can take 

iy{d{x,g)) = f{x) dxTi{dg) , (3.10) 

and denoting ^{u) := number of calls using point u, 

M{{x,gm = m + h^+9}){^) < Cyu e M}. (3.11) 

Discrete processes Free birth-and-death processes with a countable family 
of individuals G are simply the product of independent birth death processes 
labeled by each 7 e G, with birth rates ^(7) and death rate equal to the 
number of alive individuals. Such a process exists without any requirement on 
the weights 10(7); it is ergodic and its invariant distribution is the product of 
Poisson laws with mean w(7). 

The (discrete) loss networks, the contour model and the animal version of 
the random cluster model of Section 2 are processes of this form where, in 
fact, the matrix M takes only two values, and 1. That is, the interaction 
imposes a deterministic constraint. In particular, the interaction terms of the 
Ising-model, random cluster and loss networks have a simple product form 

M{l\i) - n [1-^7,^)]- (3.12) 
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Indeed, the Radon-Nikodim derivative is one for allowed configurations, hence 
the denominator in (3.3) is one. 

3.2 Graphical construction 

We proceed to the construction of the probability space where both the free 
process and the interacting birth death process will be constructed. Consider 
the countable family of random quartets {(Fj, Tj, jS'j, Zj) : i £ J^}, with Fj G 
G, Tj, Si e M, e [0, 1] such that: 

• The process {(Fj,Tj) , i e J} is a Poisson process on G x R with mean 
measure v x £; i is the Lebesgue measure in R. This process determines the 
times and type of attempted births of individuals. 

• Si is exponentially distributed with mean 1. This variable will determine 
the lifetime of the ith attempted birth. 

• Zi is uniformly distributed in (0, 1). This variable is called the flag or mark 
of the ith attempted birth and will be used together with the function M(-|-) 
to decide if the attempted birth is actually a real birth. 

Each triplet (Fj, Tj, Si) can be visualized as a cylinder of (space) basis Fj, birth 
time Ti and hfetime Si. The random set of marked cyhnders is called 



For a generic marked cylinder C = (Fx [t,t + s], z) G C, denote Birth (C) = t, 
Death (C) = t + s. Life (C) = [t,t + s], Basis (C) = F and Flag (C) = z. 

The free process The construction of C is time-translation invariant. Call 



the set of individuals forming the sections of C at time t. All attempted births 
are actual births in this case, ^til) will be at most 1 in the continuous case, 
but could be bigger in the discrete case, (^t : t G M) constitutes a stationary 



C^{(TiX[Ti,T, + S,],Zi) ,tej}. 



(3.13) 



^,(7) := #{C G C : Basis (C) = 7 ; Life (C) 3 t} 



(3.14) 
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free birth-death process with generator (3.1). The marginal law of is /jP, the 
Poisson process with intensity u. 

Likewise, one can define the free process on R*^ x [0, oo) with initial configu- 
ration of individuals 

Co := Hnl---}- (3.15) 

For this associate cylinders to the initial configuration: 

Co(eo) := {(7° X [0, + S^] , Z^) , i e (3.16) 

where Sf and are independent and independcnt-of-everything random vari- 
ables whose distributions are, respectively, Exp(l) and U{0,1). Define the 
subset of cyhnders born between and t: 

C[o,t]:= {CeC : Birth (C) e [0, t]} 

Then, the process defined at time t by 

6(7) := #{C e C[o,t] U Co(Co) : Basis (C) = 7 ; Life (C) 3 t} (3.17) 

has initial configuration ^0 and generator (3.1). 



Interacting processes The absolute continuity with respect to the free 
process, embodied in the generator (3.4), suggests a simple alteration to the 
previous construction to pass to an interacting birth-and-death process: The 
attempted births become actual births only if an additional (generally stochas- 
tic) test is passed. This test is determined by the factor M of the rate densities. 
The interacting process is, therefore, obtained as a "thinning" or "trimming" 
of the free process. 

The formalization of this intuitively simple idea is easy for finite windows, but 
more delicate for the infinite- volume process. We discuss the former case first. 
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3.2.1 Finite-volume construction 

To construct a birth-and-death process with rate density i>{dj) M(7|^), for 
individuals within a finite space-region A and for a finite time interval [to, tfin], 
one proceeds as follows: 

(1) Run the free process with rate density u starting from the initial cylinders 
Cq. If M(7|^) is deterministic — for instance forbidding individuals to 
overlap — the initial configuration is assumed to satisfy the corresponding 
constraint. 

(2) Each death happening before reaching an event of the free process causes 
the corresponding updating of ^t, by taking the corresponding individual 
out of 

(3) When the free process yields a first event (71, ti, Si, Zi), this event is 
considered an attempted birth. To decide, one looks to the set ^ti- of 
alive individuals {^ti- is equal to minus the initial individuals with 
lifetime smaller than ti). If 

zi < M(7i|et,-) (3.18) 

the cylinder is allowed to be born and the individual 71 is included in the 

configuration {j^; otherwise it is ignored and ^ti is set equal to $,ti-- 

(4) Now iterate the procedure, that is, repeat the previous two steps shifting 
subscripts 1 — > 2 and — > 1. Continue in this way until reaching an 
attempted birth beyond tfin- 



3.2.2 Two-sweep finite-volume construction 

The visualization in terms of cylinders suggests an alternative implementation 
as a two-sweep scheme: In the first sweep one generates free cylinders by 
running the free process from to to tfin, while in the second sweep a decision 
is made on which cylinders are kept and which are erased. The set of kept 
cylinders includes, by definition, all initial cylinders while successive additions 
must pass the test (3.18). We call K[o,f](A, ,^0) the resultant set of kept cylinders 
in the construction of Section 3.2.1. The configuration of the process at time 
t with initial configuration ^0 is then given by the projection of the bases of 
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the alive kept cylinders at that time: 

77^(7) = #|C G K[o,t](A, {0) with basis 7 and alive at t| . (3.19) 

3.2.3 Finite-volume time-stationary construction 

The construction can be also performed in a stationary manner for i e R. 
Indeed, since in a finite window the number of alive individuals is finite (with 
probability one), there exist random times {tj G M : j G Z}, such that (a) 
Tj ±00 for j ±00 and (b) (U.j[Tj,Tj + S.i\) fl (Uj{rj}) = 0. In words, at 
each Tj no cylinder is alive. The above selection of kept cylinders can then be 
performed independently in each of the random intervals [Ti,Tj+i). This sta- 
tionary construction is particularly useful to study properties of the invariant 
measure In fact, calling K(A) the (time stationary random) set of kept 
cylinders, the law of 

,^^(7) = #{c G K(A) with basis 7 and alive at (3.20) 
is exactly //a- 

3.2.4 Infinite-volume construction 

None of the finite-volume procedures discussed above can be directly imple- 
mented to construct the process in infinite volume. On the one hand, the 
scheme proposed for finite time intervals is not applicable to infinite vohime 
because it is not possible to decide which is the first mark in time. On the 
other hand, the stationary construction is also not feasible because in infinite 
volume there are cylinders alive at all times. This last objection, however, may 
play no role if one only focuses on a family of cylinders intersecting a finite 
set and tries to decide which of them should be erased and which ones kept. 

According to the previous discussion, to decide whether a cylinder C G C is 

kept, one has to look at the set of cylinders C (born before C and) alive at the 
birth-time of C whose basis are incompatible with the basis of C in the sense 
of (3.5). Let us call this set the first generation of ancestors of C and denote 
it Af . Once we determine which of these ancestors are ahve, the decision on 
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whether to keep C or not requires only a single application of the test (3.18). 
However, to decide which of these ancestors are alive we have to work with 
the second generation of ancestors of C, that is, with the ancestors of the 
ancestors. Recursively, we find ourselves having to deal with all generations of 
ancestors of C. Let us call the union of all generations of ancestors of C the 
clan of ancestors of C, and denote it A*" — Un>iA^, where is the set of 
ancestors in the nth generation. These sets may contain cylinders in Co('Co)- 
The procedure for deciding whether to keep or to erase C can be univocally 
defined if the clan of ancestors of C is finite. 

This picture makes it apparent that an infinite-volume process given initial 
starting conditions (i.e. for a finite time-interval) exists as long as there are no 
explosions, that is, as long as no cylinder can develop infinitely many ancestors 
in a finite time. Furthermore, there exists a unique stationary process (for 
infinite time-intervals) if all clans of ancestors are finite with (free-process) 
probability one. 

Theorem 1 

(i) If with probability one A'-' r\C[o^t] is finite for every cylinder C alive at timet, 
for any t > 0, then the birth- and- death process with the generator (3.2) and 
initial condition rjo is obtained by performing the two-sweep construction of 
Section (3.2.2) on each set (A'^n C[o,t]) U Co(r7o) o-nd taking the projections 

Vt{l) — ^jC" e K[o,t](7yo) with basis 7 and alive at i| . (3-21) 

where K.[o^t]{Vo) is the resulting set of kept cylinders, 
(a) If with probability one AP is finite for every cylinder C, then the stationary 
birth- and- death process with the generator (3.2) can he constructed for t G M 
by performing the two-sweep construction of Section (3.2.2) on each set AP 
and taking the projections 

Vtil) — T^jC* e K with basis 7 and alive at t| . (3.22) 

Moreover, the marginal distribution of rjt is the stationary measure /i. 

A proof of this theorem is presented in the Appendix. 
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4 Oriented percolation and branching processes 



To determine the conditions allowing the construction of Theorem 1, we point 
out that the relation "being ancestor of" gives rise to a model of oriented- 
percolation. Wc call it backwards oriented percolation to emphasize the fact 
that it is defined by only looking into the past. The finite-time construction 
is possible if there is no cluster with infinitely many members in a finite time 
slice, while the feasibility of the infinite-time construction requires the absence 
of a percolation cluster reaching to time — oo. As usual in oriented percolation 
problems, it is useful to work with a majorizing multitype branching process. In 
this process the offspring distribution of a cylinder C has the same (marginal) 
law as the distribution of Af , but the branches behave independently. The 
problem is then reduced to determine conditions guaranteeing the finiteness of 
the clan of branching ancestors. Hence it is sufficient to show in the finite-time 
case that the branching process docs not explode, while in the infinite-time 
case wc need to prove that the branching process is sub-critical. Factorization 
makes these tasks easier. 

Let us give sufficient conditions on the dominating branching for the different 
processes listed in Section 2. 



Discrete processes If the family of individuals G is countable, the free 
birth-and-death process is the product over 7 e G of independent marked 
Poisson processes. The construction of the interacting processes is an obvious 
adaptation of the procedure of the continuous case. For the infinite-volume 
process one relies on the properties of the backwards oriented percolation 
model of cylinders defined by the oriented bonds C — > C" if C is an ancestor 
of C, that is if the basis of C and C intersect and C is ahve when C is born. 
Let m(7, 6) be the mean number of cylinders of basis 9 in the first generation 
of a cylinder of basis 7. These are cylinders born at negative times —t and 
have a lifetime at least t, so they survive to intersect the grain born at time 
zero. Its average number is, therefore. 







00 




dt / dse 



—s 




(4.1) 



—00 
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Define m"(7, 6) as the mean number of cylinders of basis 6 incompatible with a 
cylinder of basis 7 in the n-th generation of ancestors, m" is the matrix-product 
of m by itself n times. The condition for absence of oriented percolation is 

^^m'^(7,^)<oo (4.2) 

n>l e 

for all 7. For any function "size" q : G ^ M"*", such that inf-yg(7) > 1, as in 
Lemma 5.15 of Fernandez, Ferrari and Garcia (2001), calling 

ag sup-^^g(^)m(7,^), (4.3) 



we have 

^m"(7,^) < <g(7)- (4.4) 
e 

The form of measuring this "size" depends on the process in question, but 
usually there is an obvious prescription. For instance, for the loss networks, the 
Peierls contours and the random clusters model this measure is just the length 
of the call, the number of plaquettes of the contour or the number of points 
of the cluster. The (infinite-volume) birth- and- death process — and hence the 
corresponding loss network — exists for finite time-intervals if aq < 00, while 
Qjg < 1 is a sufficient condition to be an ergodic infinite-time process. 



Area interaction point process Here m{x, .) is a measure on G = W^; 
m{x, dy) represents the rate at which cylinders with basis centered at y appear. 
Consider a germ xq, which, by space-time- invariance, can be placed at the 
origin and assumed to be born at time zero. Its ancestors are all cylinders 
whose bases involve germs located in dG :— Therefore, 
as in (4.1), 

m{x, G) = v{dG) ■ 1 . (4.5) 
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We conclude that the corresponding birth-and-death process exists for finite 
times as long as u{dG) < oo and, if furthermore, 

u{dG) < 1 (4.6) 



then there is an stationary ergodic process, absolutely continuous respect to 
the free process, having as invariant measure the corresponding point process 
of Section 2. 

The argument also works if G G G is a random set chosen independently 
of everything as in (2.3). Recall f{x) is the intensity of germs and tTx is the 
distribution of the grain centered in x. Let 

m{G, dH) = J dxf{x) I{G, H) Tr^idH) (4.7) 



the rate at which individuals H having an influence in the birth-rate of G 
appear (see (3.5) for the definition of J). This implies that the mean number 
of individuals in the first generation of ancestors of G is m{G, G). Let the 
"matrix product" be defined inductively by = m and 

m"(G', dH) := J m"-\G, dK) m{K, dH) n>l . (4.8) 

G 



As for n = 1, 'm'^{G, G) is the mean number of individuals in the n-th gener- 
ation of ancestors of G. The ergodicity of the process is imphed by 

^m"(G', G)<oo (4.9) 

n>l 



for all individual G. A sufficient condition for (4.9) is 

1 

G q{G) 



ag := sup I fix) dx J 7r,{dH) q{H) I{G, H) < 1 . (4.10) 



G 
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for some function q : B{M.'^) M"*" satisfying infGg(G') > 1. Indeed, it can be 
proven as in Lemma 5.15 in Fernandez, Ferrari and Garcia (2001) 



Strictly speaking the above statements have been rigorously proven only for 
the discretized version of the models. 



Loss networks The calls of the loss networks can be interpreted as germ- 
grains. For instance, in the one-dimensional case, the germs are the leftmost 
points of the calls and the grains are segments with random lengths. 

A particular case where one can explicitly compute the sufficient condition 
is the one-dimensional continuous loss networks of Section 3.1. Assume, in 
general, that the leftmost points of calls appear with rate f{x) and that call 
lengths are given by a distribution tt independent of x. We only require the 
latter to have a finite mean pi. Consider a germ sitting at the origin, that is 
a call stretching from the origin to the right, born at time zero. Its ancestors 
correspond to cylinders with sufficient lifetime and with bases given by either 
calls starting at negative sites and passing through the origin, or calls of arbi- 
trary length originating within the sites occupied by the initial call. Therefore, 
the aq in (4.10) for the case q{L) = 1 turns to be: 



In the homogeneous case {f{x) = k) this gives the following condition for 
ergodicity: 



A simple computation shows that choosing q{L) = max(L, 1) gives ag < 
k{p2 + pi + l), where pi and p2 are the first and second moment of the distribu- 



m'^(G,G) < q{G)a'^. 



(4.11) 




(4.12) 



K (pi -I- supL) < 1. 



(4.13) 



L 
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tion TT respectively. This gives the following sufficient condition for ergodicity 

K{p2 + Pl + l) < 1. (4.14) 



Marie (2002) improved this bound to 

«(v^ + Pi)<l- (4.15) 

We remark that to obtain these conditions it was important to consider only 
oriented percolation. The analogous conditions obtained by considering un- 
oriented percolation of cylinders are far more restrictive. 



5 Perfect simulation of invariant measures of birth-and-death pro- 
cesses 

The main issue of this section is a construction of the set A^'° formed by the 
cylinders with bases intersecting the space-time set A x {0} ("cylinders alive 
at time 0" ) and their clans of ancestors. This is a problem of simulation of 
cylinders generated by the free process. Once these clans are perfectly sim- 
ulated, it is only necessary to apply the deterministic "cleaning procedure", 
based on the test (3.18), to obtain a perfect sample of the interacting process. 
The scheme is feasible if these clans are finite with probability one, a fact valid 
under conditions like (4.6), (4.9) or ag < 1, where is defined in (4.3) for 
the discrete case and (4.10) for the area interaction process. 

We propose a non-homogeneous time-backwards construction of the clan based 
on a result proven in Section 4.5.1 of Fernandez, Ferrari and Garcia (2001). 
It is shown there that the clan of ancestors of a family of cylinders can be 
obtained combing back in time and generating births of ancestors with an 
appropriate rate. Alternatively, one could use the fact that the law of C is 
time-refiexion invariant, to generate deaths of ancestors. This is simple in the 
area-interaction process with a fixed grain, but it is not Markovian and more 
involved in the infinite case. This approach was proposed by one of the referees 
and developed by Garcia (2000). 
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For concreteness, let us discuss our scheme for individuals living in M.'^ or Z*^. 
The birth-rate of a new cylinder to be added to the clan is equal to the rate 
density of the free process multiplied by an exponential time factor ensuring 
that the ancestor has a lifespan large enough to actually be an ancestor. This 
time factor involves the time-distance to the birth of existing cylinders, which 
can be expressed through the following function. For a finite region A and a 
finite set of cylinders H, let the set of bases of the potential ancestors of H 
and A x {0} be 

G(H, A) := e G : /(Basis (C), 9) = 1, for some C e h} 

U{^eG:enA^0} (5.1) 

and for a given individual 9 e G(H, A), 

TI(H, A, 9) = min{Birth (C) : C" e H, /(Basis (C), 9) = l} (5.2) 



with the convention min0 — 0. By definition, TI(H, A, 9) < 0. 

Theorem 2 The clan A'^'^ is the limit as t ^ oo of a process A^, defined by 
the initial condition Aq = and the evolution equation 



E{^^\K,0<s<t) = 



dt 

oo 

J u{d9) J ds e-' [F{At U {9, -t, s)) - F(A,)] . (5.3) 

G(At,A) t+TI(At,A,e) 

Here F is an arbitrary function depending on a finite number of individuals 
intersecting A and we have denoted {9, —t, s) the cylinder of base 9, bom at 
time —t and with lifetime s. 

For completeness, a proof of this theorem is presented in the Appendix. For 
the free discrete loss-network (contours, random cluster) processes, 

/ u{d9)F{9) = E w{9)F{9) (5.4) 
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while for free birth-and-death processes on IR*^, 

J u{de)F{e) = j f{x)dx j ',^{dg)l{x®geG{Kt,K)]F{x,g). 

G(At,A) Rd G:, 

Notice that is a monotone process (A^ C A^^^) in which at time t only 
cylinders in G(At, A) can be included. The inclusion of a cylinder born at time 
—t requires that either (a) its basis is incompatible with that of some cylinder 
born later and its lifespan reaches the birth-time of such cylinder, or (b) its 
basis is compatible with those of all cylinders born later, but it intersects A 
and the cylinder survives up to time equal zero. The last condition is ensured 
via the convention min = in the definition of TI. 



Algorithm to construct the backwards clan of a finite region The 

combination of (5.1)/(5.3) can be translated into the following explicit al- 
gorithm. We do it first for the case of countable number of individuals and 
indicate at the end of this section how to proceed in the continuous case. To 
generate A^''^: 

(1) Start with r = and H = 0. 

(2) Let H be the current set of cylinders and r the current — min{Birth (C) : 
C e H}. For each 7 e G(H,A) generate an independent realization of 
the first time Ti(7) of the non- homogeneous Poisson process in M with 
intensity 

X^{ds) := w(7) e-^+TW.T) > t] ds . (5.5) 

Notice that ti(7) may be infinity. 

(3) Order the set {ti(7) : 7 G G(H,A)}. Let f be the infimum of this set. 
This is well defined and strictly positive because the condition a < 1 
guarantees that the total rate Z]7eG(H,A) /r+ A^((is) < 00. 

(4) If f < 00, call 7i the basis corresponding to the minimum (i.e. ti(7i) = r). 
Update r ^ f and H ^ H U {(71, -r, r + TI(H, A, 7) + where i?i 
is an exponential random variable with rate 1 independent of everything. 
In the sequel ignore the set {ti(7) : 7(7,71) = 1} (we can reuse the 
remaining ri) and go to (2). 
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(5) If r = cx) set A'^'^ = H and stop. By Theorem 2 the distribution of the 
set A^'° so generated is exactly that of the free birth-and-death process. 

If Tj are the successive times of jump of A^, then Hj = A^-. have the same 
distribution as the i-th iterate of the above algorithm. 

In the continuous case, time and space cannot be in general separated. Instead 
of steps (2) and (3) above we must consider the events (7, s) of a Poisson 
process on G x R"*" with intensity 

A((i(7, s)) = u{d-i) e-^+Ti(H,A,7) > ^1 g G(H, A)} ds . (5.6) 

For a finite window A the total rate is finite, hence these events can be well 
ordered by looking to the time coordinate. If the set of these events is not 
empty, we take f to be the minimal time coordinate (it is strictly positive 
with probability one) and denote 71 the associated individual. If the Poisson 
process with rate density (5.6) yields no event we take f — 00. We then 
continue as in (4). 

This algorithm plus the subsequent "cleaning algorithm" constitutes our per- 
fect simulation scheme. 

The cleaning algorithm Let A^'° be the clan of the cylinders whose life 
contains time and basis intersects A. The following algorithm shows how to 
construct inductively the set K^'° of kept cylinders. 

(1) Start with H = A^'° and K = (H is formed by the cyhnders to be 
tested and K by those ah^cady kept). 

(2) If H is empty go to 5. If not, order the cylinders of H by time of birth. 
Let Ci be the first of those cylinders; call 71 its basis and ri its birth- 
time. Let ^1 be the set of bases of the cyhnders in K alive at ri which are 
incompatible with the basis of Ci. Let Zi be a random variable uniformly 
distributed in [0, 1] independent of everything. 

(3) If Zi < M(7i|{i), then update: H ^ H \ {Ci}, K ^ KU{(:7i}. Go to 2. 

(4) If Zi > M(7i|^i), then update: H ^ H \ {Ci}. Go to 2. 

(5) Set K^'° = K and stop. By Theorem 1 (ii) the distribution of this clan 
K^'° is exactly that of the interacting birth-and-death process. 
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Algorithm to simulate a finite window of jj, This is the easiest part. 
Once the set K'^'^ of kept cylinders has been determined, take the configuration 
Tj defined by 

7y(^) — ^ \{C has basis 7 and hfe containing 0} (5.7) 

for 7 intersecting A. This configuration has the marginal distribution of the 
infinite- volume measure ^ on the (not necessarily finite) set Ga — pj{^^ct;enA7^0}_ 
This fact is guaranteed by Theorem 1 (ii). 



6 Errors in Perfect simulation? 

Even in finite volume, perfect simulation algorithms are subjected to error. In 
general terms, a perfect simulation algorithm of a measure on a set A" is a 
function $ : [0, 1]^ — A", such that, if {Un)nm is a sequence of i.i.d. uniform 
in [0, 1] random variables, there exists a stopping time T for ([/„) such that $ 
depends only on the first T coordinates of (C/i, U2, ■ ■ ■) and 

P($(t/i,...,t/T) e^) =M^) ■ (6-1) 

The CFTP algorithm, for instance, stops when a random value t is found such 
that the different copies of the algorithm coupled from time —t started with all 
possible initial conditions attain the same configuration at time 0. Finding t 
requires the use of a random number T{t) of uniform random variables, which 
must be less than S ="the maximum time left in order to have the results 
ready for the next congress", for instance. Thus, one actually samples from 
the distribution defined by 

p($([/i,...,c/t) e A|r<5) (6.2) 

which is different from, though as 5* — > 00 converges to, (6.1). This is the so- 
called impatient-user bias. The CFTP algorithm also permits the construction 
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of a joint realization with marginals (6.1) and (6.2) such that T < S 

implies 77 = ^. In fact, as pointed in Proposition 6.2 of Fill (1998) 



sup 

A 



p($([/i, ...,UT)eA I T <S)- i^{A) 



< 



F[T > S] 
1 - P[T > S] 



(6.3) 



In our algorithm T is determined by the number of uniform random variables 
necessary to construct the clan of the observed region A. 

When the possible sizes of the individuals 7 form an unbounded set, for in- 
stance for the Peierls contours of the Ising model, practical limitations prevent 
the inclusion of all possible sizes in the simulation. In fact the mere enumera- 
tion of the possible contours is beyond reach when more than a few dozens of 
links are involved. This is tantamount to a "space impatient-user bias": the 
user is forced to do a space cut-off that produces a bias, even when the actual 
probability for a cut event to take place is tiny. In mathematical terms, one 
actually samples from the conditioned measure 

P($(t/i, ...,UT)eA I {K <k}n{T < S}) (6.4) 



where K — "maximum perimeter of bases of cylinders in the clan" (A; = 30, 
for instance). In fact, our approach also admits a joint realization with 
marginal distributions (6.1) and (6.4) such that rj = ^ ii K < k and T < S, 
and such that F{{K > k} n {T > S}) goes to zero exponentially fast in S and 
in the cutoff of the length of the contours (30 in our example). Shghtly more 
precisely, a bound hke (6.3) holds with 

P({i^ > A;} n {T > S}) < 0(a^ X sup7r^(X > k)) . (6.5) 



This follows from the subcriticality of the majorizing branching process. For 
the Ising model, for instance, nx{K > k) — 0{e~'^^). 
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7 Conclusion 



Our algorithm offers an approadi to perfect simulations of processes with 
infinite state space. The fact that there is no coupling between different initial 
conditions, makes it a flexible tool for processes with a large state space. No 
"sandwiching processes" need to be followed; the free process is a natural 
"dominating process" in our setting. In addition, our algorithm is backed by 
a rather detailed theory that allows the estimation of various properties of 
the resulting measure, as well as possible errors. In particular, our approach 
is not free from the "impatient- user bias" , but the resulting error is relatively 
straightforward to control. 

A noteworthy feature of our approach is that the perfect simulation stage 
applies, in fact, to the free process. Interacting processes are then obtained by a 
deterministic "cleaning" . As a consequence our scheme allows the simultaneous 
simulation of all processes absolutely continuous with respect to the same free 
process. This coupled construction could be potentially useful, for instance to 
establish comparison criteria. 

The algorithm admits a further generalization more or less immediate that 
has not been pursued here: it can be applied to processes with variable death 
rate that, however, must be uniformly bounded from below by 1, say. The 
dependences in the birth and death-rates induce definitions of incompatibility 
and respective parameters a (cf. (4.10), (4.3)). A construction analogous to the 
one described in this paper can be performed but with a thinning algorithm 
that takes also into account the variable death-rates. 

In this work, the advantages of the approach have been exploited only at a 
theoretical level, where it has led to a new treatment of systems with exclu- 
sions and to better estimates of regions of existence of a number of processes. 
Berthelsen and M0ller (2001) compared it to the dominated CFTP introduced 
by Kendall and M0ller (2000). Based on simulation results, the authors show 
that the dominated CFTP is better than the algorithm based on the clan of an- 
cestors in the particular case of a Strauss process (see Equation (2.6)) defined 
on a unit square with = 100 and = (the so-called hard-core process), 
0.5 and 1 (a Poisson processes with rate 100). This is obviously the case from 
the description of the processes since the backward construction of our algo- 
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rithm stops when the dominated Poisson process regenerates and usually the 
coupling of CFTP is achieved before it in the finite case. However, it should 
be noticed that the algorithm based on the clan of ancestors was designed for 
sampling the infinite- volume process viewed in a finite window. This seems to 
be a much more interesting and challenging problem which has been studied 
by Marie (2002) for the specific case of one- dimensional loss networks with 
bounded calls. No comparison was made to other perfect simulation schemes. 

Finally, we hope that a suitable combination of our ideas with some rejection 
sampling scheme could yield a version free of the user-impatience bias. 



A Proof of Theorems 

Proof of Theorem 1 We need to show that has generator (3.4). Denote 
rit — rjt and K[0,t] the set of kept cylinders born at time zero or after time 
zero (this includes the cylinders induced by the initial configuration) and for 
F a function depending on individuals intersecting a region with finite total 
rate, write 

= J2 l{Birth (C) e[t,t + h]}[F{rjt + 5Basis(C)) - F{rjt)] 
ceK[o,t+/i] 

+ l{Ufe{C)3t,me{C)^t + h}[F{r]t-5B^,,^c))-F{r]t)] 

ceK[o,t] 

+ {other things}, (A.l) 

where {other things} refer to events with more than one Poisson mark in the 
time interval [t, t + h] for the contours in the (finite) support of F. Since the 
total rate of the Poisson marks in this set is finite, the event {other things} 
has a probability of order (/im(Supp (F), G))^, where m is defined in (4.7). 
We have 



^ l{Birth (C) e[t,t + h]} 1{C e K[0, t + h]} [F{r)t + (5Ba.is(C)) - F{r)t)] 
cec 
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= l{Birth (C) e[t,t + h]} l{Flag (C) < M(Basis (C)|77t)} 

) - nvt)] (A.2) 



To compute the second term of (A.l), observe that Life (C) is independent of 
Birth (C) and both the event {C G K[0, t]} and rit are jF^-measurable. Here JF^ 
is the cr-algebra generated by the births and deaths occurred before t. Hence 



P(Life (C) 3 t, Life (C)^t + h\ 

= P(Life {C)^t + h\ Life (C) 9 l{Life (C) 9 i} 



and 



1{C e K[0, i]} l{Life (C) 9 Life (C) ^ t + /i} 

c 

= E[^p(Life {C)^t + h\ Life (C) 9 t) 1{C e K[0, i], Life (C) 9 



X [F{r]t - 5Basis(C)) - ^(^/t)] 

Since Life(C) is exponentially distributed with mean 1, 



(A.3) 



P(Life {C)^t + h\ Life {C)3t)=h + o{h) . (A.4) 
Taking the expectation of (A.l) and substituting (A.2)-(A.4) we get 



E[F(ry,+,) - F(r/,)] 

j v{d^) e(m(7|77,) [F(77, + 5,) - F(77,)]) + o{h) 

G 

+ h Yl ^{vt{7)[F{vt-S^)-F{rit)])+o{h) (A.5) 

f:r]t{ga)>0 
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which dividing by h and taking hmit gives 



dt 



AEF{rj^). 



□ 



(A.6) 



Proof of Theorem 2 Define 



At = {C e A^'° : > Birth (C) > -t} 



A^'°nCH,0], (A.7) 



that is, the set of cyhnders in A^''^ with birth-time posterior to —t. It suffices 
to prove that the process so defined satisfies the evolution equation (5.3). 

The inclusion of a new cylinder in the time interval [t,t + h] depends on 
the existence of a birth Poisson mark in [—t — h, —t] whose corresponding 
cylinder is incompatible with some C G Aj. That is, if C is a cylinder with 
/(Basis (C), Basis (C)) = 1 for some C e At, 



The remainder o{h) is the correction related to the probability that C is not 
the only relevant cylinder born in [—t — h, —t] . Hence 



where m is defined in (4.7) and | At| stands for J2c'eAt Q'(Basis (C)), where q is 
the measure used to define a in (4.10). Since the birth-time is independent of 
the lifetime which is exponentially distributed with rate one. 



F[At+h^AuC\At^A,At>^At>,t'e[0,t)) = P{C e C : 
Birth (C) e [-t - h, -t] , Death (C) > t - TI(A, A, Basis (C))} + o{h) 



01 



(h) <{hJ2 ^i(Basis (C), G)) < h'^ \At\^ o? 




'At+h = A U C I At = A, At. = A*., e [0, t)) 
P{C e C : Birth (C) e [-t - h, -t]] 

xP(Life (C)>t- TI(A, A, Basis (C))) + o{h) 
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= /i/(Basis(C))e-*+^^(^'^'^^^^(^» + o(h) 



(A.8) 



This implies that when the configuration at time t is A, a new cyhnder with 
basis 7 is included in Aj(T) at rate 

/(7) e-*+T^(^'^'^). (A.9) 

From (A.8), as in the computation of the forward Kolmogorov equations, we 
get (5.3). This equation characterizes the law of the process At(T) as a non- 
homogeneous Markov process. □ 
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